// hnummatrix.h
#pragma once

#ifndef __HNUMMATRIX_H__
#define __HNUMMATRIX_H__

#include "hnumdef.h"

namespace harlinn
{
    namespace numerics
    {
        enum class MatrixStorageType
        {
            // column-wise, no supernode
            NC,    
            // column-wise, column-permuted, no supernode 
            // (The consecutive columns of nonzeros, after permutation, may not be stored  contiguously.) 
            NCP,   
            // row-wize, no supernode
            NR,    
            // column-wise, supernode
            SC,    
            // supernode, column-wise, permuted
            SCP,   
            // row-wise, supernode
            SR,    
            // Fortran style column-wise storage for dense matrix
            DN,     
            // distributed compressed row format
            NR_loc  
        };


        enum class MatrixDataType
        {
            // float
            S, 
            // double
            D,
            // std::complex<float>
            C,
            // std::complex<double>
            Z
        };

        




        enum class MatrixType
        {
            // general
            GE, 
            // lower triangular, unit diagonal
            TRLU,  
            // upper triangular, unit diagonal 
            TRUU,  
            // lower triangular
            TRL,   
            // upper triangular
            TRU,   
            // symmetric, store lower half 
            SYL,   
            // symmetric, store upper half 
            SYU,   
            // Hermitian, store lower half
            HEL,   
            // Hermitian, store upper half
            HEU    
        };

        /* Stype == SLU_NC (Also known as Harwell-Boeing sparse matrix format) */
        template<typename T, typename int_t = int>
        struct NCFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::NC;
            }
            // number of nonzeros in the matrix
            int_t  nnz;
            // pointer to array of nonzero values, packed by column 
            T *nzval;

            // pointer to array of row indices of the nonzeros 
            int_t  *rowind; 
            // pointer to array of beginning of columns in nzval[] and rowind[]
            int_t  *colptr; 
            // Note: 
            //  Zero-based indexing is used; colptr[] has ncol+1 entries, 
            //  the last one pointing beyond the last column, 
            //  so that colptr[ncol] = nnz. 


        };

        /* Stype == SLU_NR */
        template<typename T, typename int_t = int>
        struct NRFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::NR;
            }

            // number of nonzeros in the matrix 
            int_t  nnz;
            // pointer to array of nonzero values, packed by raw 
            T* nzval;

            // pointer to array of columns indices of the nonzeros 
            int_t  *colind;
            // pointer to array of beginning of rows in nzval[] and colind[] 
            int_t  *rowptr; 
            // Note:
		    //  Zero-based indexing is used;
		    //  rowptr[] has nrow+1 entries, the last one pointing
		    //  beyond the last row, so that rowptr[nrow] = nnz. 
        };

        /* Stype == SLU_SC */
        template<typename T, typename int_t = int>
        struct SCFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::SC;
            }

            int_t  nsuper;     /* number of supernodes, minus 1 */

            int_t  nnz;	     /* number of nonzeros in the matrix */
            T *nzval;       /* pointer to array of nonzero values, packed by column */
            int_t *nzval_colptr;/* pointer to array of beginning of columns in nzval[] */
            int_t *rowind;     /* pointer to array of compressed row indices of 
			        rectangular supernodes */
            int_t *rowind_colptr;/* pointer to array of beginning of columns in rowind[] */
            int_t *col_to_sup;   /* col_to_sup[j] is the supernode number to which column 
			        j belongs; mapping from column to supernode number. */
            int_t *sup_to_col;   /* sup_to_col[s] points to the start of the s-th 
			        supernode; mapping from supernode number to column.
		                e.g.: col_to_sup: 0 1 2 2 3 3 3 4 4 4 4 4 4 (ncol=12)
		                        sup_to_col: 0 1 2 4 7 12           (nsuper=4) */
                                /* Note:
		                Zero-based indexing is used;
		                nzval_colptr[], rowind_colptr[], col_to_sup and
		                sup_to_col[] have ncol+1 entries, the last one
		                pointing beyond the last column.
		                For col_to_sup[], only the first ncol entries are
		                defined. For sup_to_col[], only the first nsuper+2
		                entries are defined. */
        };

        /* Stype == SLU_SCP */
        template<typename T, typename int_t = int>
        struct SCPFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::SCP;
            }

            /* number of supernodes */
            int_t  nsuper; 
            /* number of nonzeros in the matrix */
            int_t  nnz;	 
            /* pointer to array of nonzero values, packed by column */
            T* nzval;
            /* nzval_colbeg[j] points to beginning of column j in nzval[] */
            int_t  *nzval_colbeg;
            /* nzval_colend[j] points to one past the last element of column j in nzval[] */
            int_t  *nzval_colend;
            /* pointer to array of compressed row indices of rectangular supernodes */
            int_t  *rowind; 
            /* rowind_colbeg[j] points to beginning of column j in rowind[] */
            int_t *rowind_colbeg;
            /* rowind_colend[j] points to one past the last element of column j in rowind[] */
            int_t *rowind_colend;
            /* col_to_sup[j] is the supernode number to which column j belongs; mapping from column to supernode. */
            int_t *col_to_sup; 
            /* sup_to_colbeg[s] points to the start of the s-th  supernode; mapping from supernode to column.*/
            int_t *sup_to_colbeg; 

            /* sup_to_colend[s] points to one past the end of the
			   s-th supernode; mapping from supernode number to
			   column.
		       e.g.: col_to_sup: 0 1 2 2 3 3 3 4 4 4 4 4 4 (ncol=12)
		               sup_to_colbeg: 0 1 2 4 7              (nsuper=4)
			        sup_to_colend: 1 2 4 7 12
             */
            int_t *sup_to_colend; 
            /* Note:
		     Zero-based indexing is used;
		     nzval_colptr[], rowind_colptr[], col_to_sup and
		     sup_to_col[] have ncol+1 entries, the last one
		     pointing beyond the last column.         
             */
        };

        /* Stype == SLU_NCP */
        template<typename T, typename int_t = int>
        struct NCPFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::NCP;
            }

            /* number of nonzeros in the matrix */
            int_t nnz;
            /* pointer to array of nonzero values, packed by column */
            T* nzval;
            /* pointer to array of row indices of the nonzeros */
            /* Note: nzval[]/rowind[] always have the same length */
            int_t *rowind;
		    
            /* colbeg[j] points to the beginning of column j in nzval[] and rowind[]  */
            int_t *colbeg;
            /* colend[j] points to one past the last element of columnj in nzval[] and rowind[]  */
            int_t *colend;
		    /* Note:
		        Zero-based indexing is used;
		        The consecutive columns of the nonzeros may not be 
		        contiguous in storage, because the matrix has been 
		        postmultiplied by a column permutation matrix. */
        };

        /* Stype == SLU_DN */
        template<typename T, typename int_t = int>
        struct DNFormat
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::DN;
            }

            /* leading dimension */
            int_t lda; 
            /* array of size lda*ncol to represent a dense matrix */
            void* nzval; 
        };

        /* Stype == SLU_NR_loc (Distributed Compressed Row Format) */
        template<typename T, typename int_t = int>
        struct NRFormatLoc
        {
            static MatrixStorageType StorageType()
            {
                return MatrixStorageType::NR_loc;
            }

            /* number of nonzeros in the local submatrix */
            int_t nnz_loc; 
            /* number of rows local to this processor */
            int_t m_loc;
            /* global index of the first row */
            int_t fst_row;
            /* pointer to array of nonzero values, packed by row */
            T* nzval;    
            /* pointer to array of beginning of rows in nzval[] and colind[]  */
            int_t *rowptr;
            /* pointer to array of column indices of the nonzeros */
            int_t *colind;
            /* Note:
			   Zero-based indexing is used;
			   rowptr[] has n_loc + 1 entries, the last one pointing
			   beyond the last row, so that rowptr[n_loc] = nnz_loc.
             */
        };


        template<typename T, typename int_t = int>
        class Matrix
        {
            class Data : public T
            {
                long referenceCount;
                MatrixType matrixType;
                // number of rows
                int_t  nrow;
                // number of columns
	            int_t  ncol;
            public:
                Data()
                    : referenceCount(1),
                      nrow(0),
                      ncol(0)
                {
                    T* ptr = this;
                    memset(ptr,0,sizeof(T));
                }

                long AddRef()
                {
                    return InterlockedIncrement(&referenceCount);
                }

                long Release()
                {
                    long result = InterlockedDecrement(&referenceCount);
                    if(!result)
                    {
                        delete this;
                    }
                    return result;
                }

                static MatrixStorageType StorageType()
                {
                    return T::StorageType();
                }


            };

            typedef Data data_t;
            mutable data_t* data;
            data_t* Data()
            {
                if(!data)
                {
                    data = new data_t();
                }
                return data;
            }

            const data_t* Data() const
            {
                if(!data)
                {
                    data = new data_t();
                }
                return data;
            }
        public:
            Matrix()
                : data(nullptr)
            {}
            Matrix(const Matrix& other)
                : data(other.data)
            {
                if(data)
                {
                    data->AddRef();
                }
            }
            Matrix(Matrix&& other)
                : data(other.data)
            {
                other.data = nullptr;
            }
            ~Matrix( )
            {
                if(data)
                {
                    data->Release();
                }
            }

            Matrix& operator = (const Matrix& other)
            {
                if(data != other.data)
                {
                    if(data)
                    {
                        data->Release();
                    }
                    data = other.data;
                    if(data)
                    {
                        data->AddRef();
                    }
                }
                return *this;
            }
            Matrix& operator = (Matrix&& other)
            {
                if(this != &other)
                {
                    if(data)
                    {
                        data->Release();
                    }
                    data = other.data;
                    if(data)
                    {
                        data->AddRef();
                    }
                }
                return *this;
            }

            static MatrixStorageType StorageType()
            {
                return T::StorageType();
            }

        };


        typedef Matrix< NCFormat< float > > MatrixNCS;
        typedef Matrix< NCFormat< double > > MatrixNCD;
        typedef Matrix< NCFormat< std::complex<float> > > MatrixNCC;
        typedef Matrix< NCFormat< std::complex<double> > > MatrixNCZ;
        typedef Matrix< NCPFormat< float > > MatrixNCPS;
        typedef Matrix< NCPFormat< double > > MatrixNCPD;
        typedef Matrix< NCPFormat< std::complex<float> > > MatrixNCPC;
        typedef Matrix< NCPFormat< std::complex<double> > > MatrixNCPZ;
        typedef Matrix< NRFormat< float > > MatrixNRS;
        typedef Matrix< NRFormat< double > > MatrixNRD;
        typedef Matrix< NRFormat< std::complex<float> > > MatrixNRC;
        typedef Matrix< NRFormat< std::complex<double> > > MatrixNRZ;
        typedef Matrix< SCFormat< float > > MatrixSCS;
        typedef Matrix< SCFormat< double > > MatrixSCD;
        typedef Matrix< SCFormat< std::complex<float> > > MatrixSCC;
        typedef Matrix< SCFormat< std::complex<double> > > MatrixSCZ;
        typedef Matrix< SCPFormat< float > > MatrixSCPS;
        typedef Matrix< SCPFormat< double > > MatrixSCPD;
        typedef Matrix< SCPFormat< std::complex<float> > > MatrixSCPC;
        typedef Matrix< SCPFormat< std::complex<double> > > MatrixSCPZ;
        typedef Matrix< DNFormat< float > > MatrixDNS;
        typedef Matrix< DNFormat< double > > MatrixDND;
        typedef Matrix< DNFormat< std::complex<float> > > MatrixDNC;
        typedef Matrix< DNFormat< std::complex<double> > > MatrixDNZ;
        typedef Matrix< NRFormatLoc< float > > MatrixNRLocS;
        typedef Matrix< NRFormatLoc< double > > MatrixNRLocD;
        typedef Matrix< NRFormatLoc< std::complex<float> > > MatrixNRLocC;
        typedef Matrix< NRFormatLoc< std::complex<double> > > MatrixNRLocZ;


        




    };
};


#endif // __HNUMMATRIX_H__

